######################################################################################
##Replication Code for ###############################################################
##Decision by Design: Leaders, Bureaucracies, and International Crisis Performance ###
##Tyler Jost #########################################################################
##International Studies Quarterly ####################################################
##(3) Appendix Analysis ##############################################################
######################################################################################

rm(list = ls())
require(here)
require(interplot)
require(sandwich)
require(clubSandwich)
require(lmtest)
require(stargazer)
require(xtable)
require(dfidx)
require(mlogit)


#########################################################################
##Institutions: Summary Statistics ######################################
#########################################################################

load("nsi_dataset_v1.RData")
vars <- c("type", "info_up", "info_horizontal", "nsdb_total", "nsdb_staff", "nsdb_mod_rep", "nsdb_dip_rep", "nsdb_fio_rep", 
          "nscb_total", "nscb_staff", "nscb_mod_rep", "nscb_dip_rep", "nscb_fio_rep")
d.summary <- subset(d, select = vars)

type_row <- data.frame(
  Variable = "type",
  Observations = sum(!is.na(d.summary$type)),
  Missing = sum(is.na(d.summary$type)),
  Mean = NA,
  SD = NA,
  Min = NA,
  Max = NA
)

numeric_data <- d.summary[sapply(d.summary, is.numeric)]
numeric_summary <- data.frame(
  Variable = names(numeric_data),
  Observations = sapply(numeric_data, function(x) sum(!is.na(x))),
  Mean = sapply(numeric_data, function(x) mean(x, na.rm = TRUE)),
  SD = sapply(numeric_data, function(x) sd(x, na.rm = TRUE)),
  Min = sapply(numeric_data, function(x) min(x, na.rm = TRUE)),
  Max = sapply(numeric_data, function(x) max(x, na.rm = TRUE)),
  Missing = sapply(numeric_data, function(x) sum(is.na(x)))
)

summary <- rbind(type_row, numeric_summary)
summary$Variable <- c("Institutional Type",
                      "Information Search Score",
                      "Bureaucratic Access Score",
                      "NSDB: Size",
                      "NSDB: Staff",
                      "NSDB: Defense Representation",
                      "NSDB: Diplopmatic Representation",
                      "NSDB: Foreign Intelligence Representation",
                      "NSCB: Size",
                      "NSCB: Staff",
                      "NSCB: Defense Representation",
                      "NSCB: Diplopmatic Representation",
                      "NSCB: Foreign Intelligence Representation")
print(xtable(summary), include.rownames = F)

#########################################################################
##Appendix Tables - Table 1 #############################################
#########################################################################

summary_table <- rbind(type_row, numeric_summary)
summary_table[,1] <- c("Institutional Type", "Information Search Score", "Bureaucratic Access Score", "NSDB: Size", "NSDB: Staff", 
                       "NSDB: Defense Representation", "NSDB: Diplomatic Representation", "NSDB: Foreign Intelligence Representation", 
                       "NSCB: Size", "NSCB: Staff", "NSCB: Defense Representation", "NSCB: Diplomatic Representation", "NSCB: Foreign Intelligence Representation")
summary_table[, c("Mean", "SD", "Min", "Max")] <- round(summary_table[, c("Mean", "SD", "Min", "Max")], 4)
print(xtable(summary_table), include.rownames = F)

#########################################################################
##Appendix Tables - Table 2 #############################################
#########################################################################

rm(list = ls())
load("analysis.RData")
cross <- table(d$crisis_init_failA[d$crisis_initA==1], d$type[d$crisis_initA==1])
prop_cross <- round(prop.table(cross, 2), 2) * 100
combined <- matrix(
  paste(cross, " (", prop_cross, "%)", sep = ""),
  nrow = nrow(cross),
  ncol = ncol(cross)
)
print(xtable(combined), include.rownames = F)

#########################################################################
##Appendix Models - Parochial/CMR Models ################################
#########################################################################

R2logit <- function(y,model){
  R2<- 1-(model$deviance/model$null.deviance)
  return(R2)
}

crse <- function(x, df){
  d.cse <- na.omit(df[, c("crisno", all.vars(formula(x)))])
  fc <- d.cse$crisno
  m <- length(unique(fc))
  k <- length(coef(x))
  u <- estfun(x)
  u.clust <- matrix(NA, nrow=m, ncol=k) 
  for(j in 1:k){
    u.clust[,j] <- tapply(u[,j], fc, sum)
  }
  bread <-vcov(x)
  c.vcov <- bread %*% ((m / (m-1)) * t(u.clust) %*% (u.clust)) %*% bread
  result <<- coeftest(x, c.vcov)
  return(result)
}

parochial1 <- glm(crisis_init_failA ~ nsdb_dip_rep + nsdb_mod_rep + nsdb_fio_rep + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop + as.factor(statea), 
            family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
pseudoR.p1 <- round(R2logit(crisis_init_failA, parochial1),3)
obs.p1 <- nobs(parochial1)
parochial1 <- crse(parochial1, d[d$crisis_initA==1,])

parochial2 <- glm(crisis_init_failA ~ nscb_dip_rep + nscb_mod_rep + nscb_fio_rep + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop + as.factor(statea), 
            family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
pseudoR.p2 <- round(R2logit(crisis_init_failA, parochial2),3)
obs.p2 <- nobs(parochial2)
parochial2 <- crse(parochial2, d[d$crisis_initA==1,])

##Create binary variable for military representation on decision-making or coordination bodies
d$nsdb_mil_rep_b <- NA
d$nsdb_mil_rep_b <- ifelse((d$nsdb_mil_rep==0 | d$nsdb_mil_rep==-8) & (d$nscb_mil_rep==0 | d$nscb_mil_rep==-8),0,1)

cmr <- glm(crisis_init_failA ~ type + nsdb_mil_rep_b + age + gender + milservice + combat + hos_experience + gwf_personal + gwf_military + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop + as.factor(statea), 
            family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
pseudoR.cmr <- round(R2logit(crisis_init_failA, cmr),3)
obs.cmr <- nobs(cmr)
cmr <- crse(cmr, d[d$crisis_initA==1,])

stargazer(parochial1, parochial2, cmr, type = "latex", omit=c("Constant", "statea"),
          omit.stat=c("f", "ser", "aic", "ll"),
          label=c("table: parochial-cmr"),
          covariate.labels = c("Decision-Making: Diplomatic Representation", "Decision-Making: Defense Representation", "Decision-Making: Intelligence Representation", "Coordination: Diplomatic Representation", "Coordination: Defense Representation", "Coordination: Intelligence Representation", 
                               "Siloed", "Fragmented", "Dictatorial", "Military Represenation", "Leader Age", "Male", "Military Experience", "Combat Experience", "Time in Office", "Military Regime", "Personalist Regime", 
                               "Regime Transition", "State A Democraticness", "State B Democraticness", "State A Capabilities", "State B Capabilities",
                               "Relative Capabilities", "Alliance Score"),
          add.lines=list(c("Observations", obs.p1, obs.p2, obs.cmr),
                         c("McFadden Pseudo-R$^2$", pseudoR.p1, pseudoR.p2, pseudoR.cmr),
                         c("Cluster Robust SEs", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                         c("Country Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$")),
          dep.var.labels = c("Failure to Achieve Goals"),
          title = "Parochial Interests and Civil-Military Relations",
          notes.append=FALSE, no.space=T, single.row = T)

#########################################################################
##Appendix Models - Pooled models #######################################
#########################################################################

##Model 1:
m1 <- glm(crisis_init_failA ~ type, 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
pseudoR2.1 <- round(R2logit(crisis_init_failA, m1),3)
obs1 <- nobs(m1)
m1 <- crse(m1, d[d$crisis_initA==1,])

##Model 2:
m2 <- glm(crisis_init_failA ~ type + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop, 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
pseudoR2.2 <-round(R2logit(crisis_init_failA, m2),3)
obs2 <- nobs(m2)
m2 <- crse(m2, d[d$crisis_initA==1,])

stargazer(m1, m2, type = "latex", omit=c("statea"),
          omit.stat=c("f", "ser", "aic", "ll"),
          label=c("table: pooled"),
          covariate.labels = c("Siloed", "Fragmented", "Dictatorial", "Leader Age", "Male", "Military Experience", "Combat Experience", "Time in Office", "Military Regime",
                               "Personalist Regime", "Regime Transition", "State A Democraticness", "State B Democraticness", "State A Capabilities", "State B Capabilities",
                               "Relative Capabilities", "Alliance Score"),
          add.lines=list(c("Observations", obs1, obs2),
                         c("McFadden Pseudo-R$^2$", pseudoR2.1, pseudoR2.2),
                         c("Cluster Robust SEs", "$\\checkmark$", "$\\checkmark$")),
          dep.var.labels = c("Failure to Achieve Goals"),
          title = "Pooled Models",
          notes.append=FALSE, no.space=T, single.row = T)


#########################################################################
##Appendix Models - Trade and GDP #######################################
#########################################################################

##Model 1: 
econ1 <- glm(crisis_init_failA ~ type + LrealgdpA + Ltrade + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop, 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
pseudoR.econ1 <-round(R2logit(crisis_init_failA, econ1),3)
obs.econ1 <- nobs(econ1)
econ1 <- crse(econ1, d[d$crisis_initA==1,])

##Model 2:
econ2 <- glm(crisis_init_failA ~ type + LrealgdpA + Ltrade + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop + as.factor(statea), 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
pseudoR.econ2 <- round(R2logit(crisis_init_failA, econ2),3)
obs.econ2 <- nobs(econ2)
econ2 <- crse(econ2, d[d$crisis_initA==1,])

stargazer(econ1, econ2, type = "latex", omit=c("statea"),
          omit.stat=c("f", "ser", "ll", "aic"),
          label=c("table: trade_gdp"),
          covariate.labels = c("Siloed", "Fragmented", "Dictatorial", "GDP", "Trade", "Leader Age", "Male", "Military Experience", "Combat Experience", "Time in Office", "Military Regime",
                               "Personalist Regime", "Regime Transition", "State A Democraticness", "State B Democraticness", "State A Capabilities", "State B Capabilities",
                               "Relative Capabilities", "Alliance Score"),
          add.lines=list(c("Observations", obs.econ1, obs.econ2),
                         c("McFadden Pseudo-R$^2$", pseudoR.econ1, pseudoR.econ2),
                         c("Cluster Robust SEs", "$\\checkmark$", "$\\checkmark$"),
                         c("Country Fixed Effects", "\\checkmark$" , "$\\checkmark$")),
          dep.var.labels = c("Failure to Achieve Goals"),
          title = "Including Trade and GDP",
          notes.append=FALSE, no.space=T, single.row = T,
          notes = c(""))

#########################################################################
##Appendix Models - Excluding the United States #########################
#########################################################################

##Model 1: 
m1 <- glm(crisis_init_failA ~ type, 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$statea!=2,])
pseudoR2.1 <-round(R2logit(crisis_init_failA, m1),3)
obs1 <- nobs(m1)
m1 <- crse(m1, d[d$crisis_initA==1 & d$statea!=2,])

##Model 2:
m2 <- glm(crisis_init_failA ~ type + as.factor(statea), 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$statea!=2,])
pseudoR2.2 <- round(R2logit(crisis_init_failA, m2),3)
obs2 <- nobs(m2)
m2 <- crse(m2, d[d$crisis_initA==1 & d$statea!=2,])

##Model 3: 
m3 <- glm(crisis_init_failA ~ type + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop, 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$statea!=2,])
pseudoR2.3 <-round(R2logit(crisis_init_failA, m3),3)
obs3 <- nobs(m3)
m3 <- crse(m3, d[d$crisis_initA==1 & d$statea!=2,])

##Model 4:
m4 <- glm(crisis_init_failA ~ type + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop + as.factor(statea), 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$statea!=2,])
pseudoR2.4 <- round(R2logit(crisis_init_failA, m4),3)
obs4 <- nobs(m4)
m4 <- crse(m4, d[d$crisis_initA==1 & d$statea!=2,])

stargazer(m1, m2, m3, m4, type = "latex", omit=c("statea"),
          omit.stat=c("f", "ser", "aic", "ll"),
          label=c("table: drop_us"),
          covariate.labels = c("Siloed", "Fragmented", "Dictatorial", "Leader Age", "Male", "Military Experience", "Combat Experience", "Time in Office", "Military Regime",
                               "Personalist Regime", "Regime Transition", "State A Democraticness", "State B Democraticness", "State A Capabilities", "State B Capabilities",
                               "Relative Capabilities", "Alliance Score"),
          add.lines=list(c("Observations", obs1, obs2, obs3, obs4),
                         c("McFadden Pseudo-R$^2$", pseudoR2.1, pseudoR2.2, pseudoR2.3, pseudoR2.4),
                         c("Cluster Robust SEs", "$\\checkmark$", "", "$\\checkmark$"),
                         c("Country Fixed Effects", "" , "$\\checkmark$")),
          dep.var.labels = c("Failure to Achieve Goals"),
          title = "Models Excluding the United States",
          notes.append=FALSE, no.space=T, single.row = T,
          notes = c(""))

#########################################################################
##Appendix Models - Excluding the Soviet Union/Russia ###################
#########################################################################

##Model 1:
m1 <- glm(crisis_init_failA ~ type, 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$statea!=365,])
pseudoR2.1 <-round(R2logit(crisis_init_failA, m1),3)
obs1 <- nobs(m1)
m1 <- crse(m1, d[d$crisis_initA==1 & d$statea!=365,])

##Model 2:
m2 <- glm(crisis_init_failA ~ type + as.factor(statea), 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$statea!=365,])
pseudoR2.2 <- round(R2logit(crisis_init_failA, m2),3)
obs2 <- nobs(m2)
m2 <- crse(m2, d[d$crisis_initA==1 & d$statea!=365,])

##Model 3:
m3 <- glm(crisis_init_failA ~ type + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop, 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$statea!=365,])
pseudoR2.3 <-round(R2logit(crisis_init_failA, m3),3)
obs3 <- nobs(m3)
m3 <- crse(m3, d[d$crisis_initA==1 & d$statea!=365,])

##Model 4:
m4 <- glm(crisis_init_failA ~ type + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop + as.factor(statea), 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$statea!=365,])
pseudoR2.4 <- round(R2logit(crisis_init_failA, m4),3)
obs4 <- nobs(m4)
m4 <- crse(m4, d[d$crisis_initA==1 & d$statea!=365,])

stargazer(m1, m2, m3, m4, type = "latex", omit=c("statea"),
          omit.stat=c("f", "ser", "aic", "ll"),
          label=c("table: drop_ussr"),
          covariate.labels = c("Siloed", "Fragmented", "Dictatorial", "Leader Age", "Male", "Military Experience", "Combat Experience", "Time in Office", "Military Regime",
                               "Personalist Regime", "Regime Transition", "State A Democraticness", "State B Democraticness", "State A Capabilities", "State B Capabilities",
                               "Relative Capabilities", "Alliance Score"),
          add.lines=list(c("Observations", obs1, obs2, obs3, obs4),
                         c("McFadden Pseudo-R$^2$", pseudoR2.1, pseudoR2.2, pseudoR2.3, pseudoR2.4),
                         c("Cluster Robust SEs", "$\\checkmark$", "", "$\\checkmark$"),
                         c("Country Fixed Effects", "" , "$\\checkmark$")),
          dep.var.labels = c("Failure to Achieve Goals"),
          title = "Models Excluding the Soviet Union/Russia",
          notes.append=FALSE, no.space=T, single.row = T,
          notes = c(""))

#########################################################################
##Appendix Models - Excluding China #####################################
#########################################################################

##Model 1:
m1 <- glm(crisis_init_failA ~ type, 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$statea!=710,])
pseudoR2.1 <-round(R2logit(crisis_init_failA, m1),3)
obs1 <- nobs(m1)
m1 <- crse(m1, d[d$crisis_initA==1 & d$statea!=710,])

##Model 2:
m2 <- glm(crisis_init_failA ~ type + as.factor(statea), 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$statea!=710,])
pseudoR2.2 <- round(R2logit(crisis_init_failA, m2),3)
obs2 <- nobs(m2)
m2 <- crse(m2, d[d$crisis_initA==1 & d$statea!=710,])

##Model 3:
m3 <- glm(crisis_init_failA ~ type + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop, 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$statea!=710,])
pseudoR2.3 <-round(R2logit(crisis_init_failA, m3),3)
obs3 <- nobs(m3)
m3 <- crse(m3, d[d$crisis_initA==1 & d$statea!=710,])

##Model 4:
m4 <- glm(crisis_init_failA ~ type + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop + as.factor(statea), 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$statea!=710,])
pseudoR2.4 <- round(R2logit(crisis_init_failA, m4),3)
obs4 <- nobs(m4)
m3 <- crse(m4, d[d$crisis_initA==1 & d$statea!=710,])

stargazer(m1, m2, m3, m4, type = "latex", omit=c("statea"),
          omit.stat=c("f", "ser", "aic", "ll"),
          label=c("table: drop_prc"),
          covariate.labels = c("Siloed", "Fragmented", "Dictatorial", "Leader Age", "Male", "Military Experience", "Combat Experience", "Time in Office", "Military Regime",
                               "Personalist Regime", "Regime Transition", "State A Democraticness", "State B Democraticness", "State A Capabilities", "State B Capabilities",
                               "Relative Capabilities", "Alliance Score"),
          add.lines=list(c("Observations", obs1, obs2, obs3, obs4),
                         c("McFadden Pseudo-R$^2$", pseudoR2.1, pseudoR2.2, pseudoR2.3, pseudoR2.4),
                         c("Cluster Robust SEs", "$\\checkmark$", "", "$\\checkmark$"),
                         c("Country Fixed Effects", "" , "$\\checkmark$")),
          dep.var.labels = c("Failure to Achieve Goals"),
          title = "Models Excluding China",
          notes.append=FALSE, no.space=T, single.row = T,
          notes = c(""))

#########################################################################
##Appendix Models - Major Power #########################################
#########################################################################

##Model 1:
m1 <- glm(crisis_init_failA ~ type + majora + majora, 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
pseudoR2.1 <-round(R2logit(crisis_init_failA, m1),3)
obs1 <- nobs(m1)
m1 <- crse(m1, d[d$crisis_initA==1,])

##Model 2:
m2 <- glm(crisis_init_failA ~ type + majora + majora + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop, 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
pseudoR2.2 <-round(R2logit(crisis_init_failA, m2),3)
obs2 <- nobs(m2)
m2 <- crse(m2, d[d$crisis_initA==1,])

stargazer(m1, m2, type = "latex", omit=c("statea"),
          omit.stat=c("f", "ser", "ll", "aic"),
          label=c("table: major"),
          covariate.labels = c("Siloed", "Fragmented", "Dictatorial", "Major Power", "Leader Age", "Male", "Military Experience", "Combat Experience", "Time in Office", "Military Regime",
                               "Personalist Regime", "Regime Transition", "State A Democraticness", "State B Democraticness", "State A Capabilities", "State B Capabilities",
                               "Relative Capabilities", "Alliance Score"),
          add.lines=list(c("Observations", obs1, obs2),
                         c("McFadden Pseudo-R$^2$", pseudoR2.1, pseudoR2.2),
                         c("Cluster Robust SEs", "$\\checkmark$", "$\\checkmark$")),
          dep.var.labels = c("Failure to Achieve Goals"),
          title = "Models Including Control for Major Power Status",
          notes.append=FALSE, no.space=T, single.row = T,
          notes = c(""))

#########################################################################
##Appendix Models - Controlling for the Cold War ########################
#########################################################################

##Model 1: 
m1 <- glm(crisis_init_failA ~ type + cw + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop, 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
pseudoR2.1 <-round(R2logit(crisis_init_failA, m1),3)
obs1 <- nobs(m1)
m1 <- crse(m1, d[d$crisis_initA==1,])

##Model 2: 
m2 <- glm(crisis_init_failA ~ type + cw + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop + as.factor(statea), 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
pseudoR2.2 <- round(R2logit(crisis_init_failA, m2),3)
obs2 <- nobs(m2)
m2 <- crse(m2, d[d$crisis_initA==1,])

stargazer(m1, m2, type = "latex", omit=c("statea"),""
          omit.stat=c("f", "ser", "aic", "ll"),
          label=c("table: drop_us"),
          covariate.labels = c("Siloed", "Fragmented", "Dictatorial", "Cold War", "Leader Age", "Male", "Military Experience", "Combat Experience", "Time in Office", "Military Regime",
                               "Personalist Regime", "Regime Transition", "State A Democraticness", "State B Democraticness", "State A Capabilities", "State B Capabilities",
                               "Relative Capabilities", "Alliance Score"),
          add.lines=list(c("Observations", obs1, obs2),
                         c("McFadden Pseudo-R$^2$", pseudoR2.1, pseudoR2.2, pseudoR2.3, pseudoR2.4),
                         c("Cluster Robust SEs", "$\\checkmark$", "\\checkmark"),
                         c("Country Fixed Effects", "" , "$\\checkmark$")),
          dep.var.labels = c("Failure to Achieve Goals"),
          title = "Controlling for the Cold War",
          notes.append=FALSE, no.space=T, single.row = T,
          notes = c(""))

#########################################################################
##Appendix Models - Results within Democracies ##########################
#########################################################################

##Recode dictatorial as fragmented
table(d$type[d$crisis_initA==1 & d$polity2A>=7])
d$mergedtype <- as.factor(d$type)
d$mergedtype[d$mergedtype=="Dictatorial"] <- "Fragmented"

##Model 1: 
m1 <- glm(crisis_init_failA ~ mergedtype, 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$polity2A>=7,])
pseudoR2.1 <-round(R2logit(crisis_init_failA, m1),3)
obs1 <- nobs(m1)
m1 <- crse(m1, d[d$crisis_initA==1 & d$polity2A>=7,])

##Model 2: 
m2 <- glm(crisis_init_failA ~ mergedtype + as.factor(statea), 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$polity2A>=7,])
pseudoR2.2 <- round(R2logit(crisis_init_failA, m2),3)
obs2 <- nobs(m2)
m2 <- crse(m2, d[d$crisis_initA==1 & d$polity2A>=7,])

##Model 3: 
m3 <- glm(crisis_init_failA ~ mergedtype + age + gender + milservice + combat + hos_experience + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop, 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$polity2A>=7,])
pseudoR2.3 <-round(R2logit(crisis_init_failA, m3),3)
obs3 <- nobs(m3)
m3 <- crse(m3, d[d$crisis_initA==1 & d$polity2A>=7,])

##Model 4: 
m4 <- glm(crisis_init_failA ~ mergedtype + age + gender + milservice + combat + hos_experience + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop + as.factor(statea), 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$polity2A>=7,])
pseudoR2.4 <- round(R2logit(crisis_init_failA, m4),3)
obs4 <- nobs(m4)
m4 <- crse(m4, d[d$crisis_initA==1 & d$polity2A>=7,])

stargazer(m1, m2, m3, m4, type = "latex", omit=c("statea"),
          omit.stat=c("f", "ser", "aic", "ll"),
          label=c("table: democracies"),
          covariate.labels = c("Siloed", "Fragmented", "Leader Age", "Male", "Military Experience", "Combat Experience", "Time in Office", 
                               "State A Democraticness", "State B Democraticness", "State A Capabilities", "State B Capabilities",
                               "Relative Capabilities", "Alliance Score"),
          add.lines=list(c("Observations", obs1, obs2, obs3, obs4),
                         c("McFadden Pseudo-R$^2$", pseudoR2.1, pseudoR2.2, pseudoR2.3, pseudoR2.4),
                         c("Cluster Robust SEs", "$\\checkmark$", "", "$\\checkmark$", ""),
                         c("Country Fixed Effects", "" , "$\\checkmark$", "" , "$\\checkmark$")),
          dep.var.labels = c("Failure to Achieve Goals"),
          title = "Results within Democracies (Polity Score <= 7)",
          notes.append=FALSE, no.space=T, single.row = T,
          notes = c(""))

#########################################################################
##Appendix Models - Results within Non-Democracies ######################
#########################################################################

##Model 1:
m1 <- glm(crisis_init_failA ~ type, 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$polity2A<7,])
pseudoR2.1 <-round(R2logit(crisis_init_failA, m1),3)
obs1 <- nobs(m1)
m1 <- crse(m1, d[d$crisis_initA==1 & d$polity2A<7,])

##Model 2: 
m2 <- glm(crisis_init_failA ~ type + as.factor(statea), 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$polity2A<7,])
pseudoR2.2 <- round(R2logit(crisis_init_failA, m2),3)
obs2 <- nobs(m2)
m2 <- crse(m2, d[d$crisis_initA==1 & d$polity2A<7,])

##Model 3:
m3 <- glm(crisis_init_failA ~ type + age + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop, 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$polity2A<7,])
pseudoR2.3 <-round(R2logit(crisis_init_failA, m3),3)
obs3 <- nobs(m3)
m3 <- crse(m3, d[d$crisis_initA==1 & d$polity2A<7,])

##Model 4: 
m4 <- glm(crisis_init_failA ~ type + age + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop + as.factor(statea), 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$polity2A<7,])
pseudoR2.4 <- round(R2logit(crisis_init_failA, m4),3)
obs4 <- nobs(m4)
m4 <- crse(m4, d[d$crisis_initA==1 & d$polity2A<7,])

stargazer(m1, m2, m3, m4, type = "latex", omit=c("statea"),
          omit.stat=c("f", "ser", "aic", "ll"),
          label=c("table: nondemocracies"),
          covariate.labels = c("Siloed", "Fragmented", "Dictatorial", "Leader Age", "Military Experience", "Combat Experience", "Time in Office", "Military Regime",
                               "Personalist Regime", "Regime Transition", "State A Democraticness", "State B Democraticness", "State A Capabilities", "State B Capabilities",
                               "Relative Capabilities", "Alliance Score"),
          add.lines=list(c("Observations", obs1, obs2, obs3, obs4),
                         c("McFadden Pseudo-R$^2$", pseudoR2.1, pseudoR2.2, pseudoR2.3, pseudoR2.4),
                         c("Cluster Robust SEs", "$\\checkmark$", "$\\checkmark$", ""),
                         c("Country Fixed Effects", "" , "$\\checkmark$")),
          dep.var.labels = c("Failure to Achieve Goals"),
          title = "Results within Non-Democracies (Polity Score $<$ 7)",
          notes.append=FALSE, no.space=T, single.row = T,
          notes = c(""))

#########################################################################
##Appendix Models - Dropping Transitioning Regimes ######################
#########################################################################

##Model 1: 
m1 <- glm(crisis_init_failA ~ type, 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$regtransA!=1,])
pseudoR2.1 <-round(R2logit(crisis_init_failA, m1),3)
obs1 <- nobs(m1)
m1 <- crse(m1, d[d$crisis_initA==1 & d$regtransA!=1,])

##Model 2: 
m2 <- glm(crisis_init_failA ~ type + as.factor(statea), 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$regtransA!=1,])
pseudoR2.2 <- round(R2logit(crisis_init_failA, m2),3)
obs2 <- nobs(m2)
m2 <- crse(m2, d[d$crisis_initA==1 & d$regtransA!=1,])

##Model 3: 
m3 <- glm(crisis_init_failA ~ type + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop, 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$regtransA!=1,])
pseudoR2.3 <-round(R2logit(crisis_init_failA, m3),3)
obs3 <- nobs(m3)
m3 <- crse(m3, d[d$crisis_initA==1 & d$regtransA!=1,])

##Model 4: 
m4 <- glm(crisis_init_failA ~ type + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop + as.factor(statea), 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1 & d$regtransA!=1,])
pseudoR2.4 <- round(R2logit(crisis_init_failA, m4),3)
obs4 <- nobs(m4)
m4 <- crse(m4, d[d$crisis_initA==1 & d$regtransA!=1,])

stargazer(m1, m2, m3, m4, type = "latex", omit=c("statea"),
          omit.stat=c("f", "ser", "aic", "ll"),
          label=c("table: drop_transition"),
          covariate.labels = c("Siloed", "Fragmented", "Dictatorial", "Leader Age", "Male", "Military Experience", "Combat Experience", "Time in Office", "Military Regime",
                               "Personalist Regime", "State A Democraticness", "State B Democraticness", "State A Capabilities", "State B Capabilities",
                               "Relative Capabilities", "Alliance Score"),
          add.lines=list(c("Observations", obs1, obs2, obs3, obs4),
                         c("McFadden Pseudo-R$^2$", pseudoR2.1, pseudoR2.2, pseudoR2.3, pseudoR2.4),
                         c("Cluster Robust SEs", "$\\checkmark$", "$\\checkmark$"),
                         c("Country Fixed Effects", "" , "$\\checkmark$")),
          dep.var.labels = c("Failure to Achieve Goals"),
          title = "Dropping Transitioning Regimes",
          notes.append=FALSE, no.space=T, single.row = T,
          notes = c(""))


#########################################################################
##Appendix Models - Dyad Fixed Effects ##################################
#########################################################################

d$dyadid <- as.numeric(factor(paste(d$statea, "-", d$stateb, sep = "")))

##Model 2: 
m1 <- glm(crisis_init_failA ~ type + as.factor(dyadid), 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
pseudoR2.1 <- round(R2logit(crisis_init_failA, m1),3)
obs1 <- nobs(m1)
m1 <- crse(m1, d[d$crisis_initA==1,])

##Model 4:
m2 <- glm(crisis_init_failA ~ type + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop + as.factor(dyadid), 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
pseudoR2.2 <- round(R2logit(crisis_init_failA, m2),3)
obs2 <- nobs(m2) #542
m2 <- crse(m2, d[d$crisis_initA==1,])

stargazer(m1, m2, type = "latex", omit=c("Constant", "dyadid"),
          omit.stat=c("f", "ser", "ll", "aic"),
          label=c("table: dyad_fe"),
          covariate.labels = c("Siloed", "Fragmented", "Dictatorial", "Leader Age", "Male", "Military Experience", "Combat Experience", "Time in Office", "Military Regime",
                               "Personalist Regime", "Regime Transition", "State A Democraticness", "State B Democraticness", "State A Capabilities", "State B Capabilities",
                               "Relative Capabilities", "Alliance Score"),
          add.lines=list(c("Observations", obs1, obs2),
                         c("McFadden Pseudo-R$^2$", pseudoR2.1, pseudoR2.2),
                         c("Cluster Robust SEs", "$\\checkmark$", "$\\checkmark$"),
                         c("Dyad Fixed Effects", "$\\checkmark$", "$\\checkmark$")),
          dep.var.labels = c("Failure to Achieve Goals"),
          title = "National Security Institutions and International Crisis Performance",
          notes.append=FALSE, no.space=T, single.row = T,
          notes = c("Robust standard errors are clustered by crisis.$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}"))

#########################################################################
##Adversary Initiated ###################################################
#########################################################################

d$crisis_noinitA <- ifelse(d$crisis_initA==0,1,0)
d$crisis_noinit_failA <- ifelse(d$crisis_noinitA==1 & d$crisis_failA==1,1,0)

##Model 1:
m1 <- glm(crisis_noinit_failA ~ type, 
          family=binomial(link = "logit"), data = d[d$crisis_noinitA==1,])
pseudoR2.1 <- round(R2logit(crisis_init_failA, m1),3)
obs1 <- nobs(m1)
m1 <- crse(m1, d[d$crisis_noinitA==1,])

##Model 2:
m2 <- glm(crisis_noinit_failA ~ type + as.factor(statea), 
          family=binomial(link = "logit"), data = d[d$crisis_noinitA==1,])
pseudoR2.2 <- round(R2logit(crisis_noinit_failA, m2),3)
obs2 <- nobs(m2) #1147
m2 <- crse(m2, d[d$crisis_noinitA==1,])

##Model 3:
m3 <- glm(crisis_noinit_failA ~ type + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop, 
          family=binomial(link = "logit"), data = d[d$crisis_noinitA==1,])
pseudoR2.3 <-round(R2logit(crisis_noinit_failA, m3),3)
obs3 <- nobs(m3)
m3 <- crse(m3, d[d$crisis_noinitA==1,])

##Model 4:
m4 <- glm(crisis_noinit_failA ~ type + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop + as.factor(statea), 
          family=binomial(link = "logit"), data = d[d$crisis_noinitA==1,])
pseudoR2.4 <- round(R2logit(crisis_noinit_failA, m4),3)
obs4 <- nobs(m4)
m4 <- crse(m4, d[d$crisis_noinitA==1,])

stargazer(m1, m2, m3, m4, type = "latex", omit=c("statea"),
          omit.stat=c("f", "ser", "ll", "aic"),
          label=c("table: adversary_initiated"),
          covariate.labels = c("Siloed", "Fragmented", "Dictatorial", "Leader Age", "Male", "Military Experience", "Combat Experience", "Time in Office", "Military Regime",
                               "Personalist Regime", "Regime Transition", "State A Democraticness", "State B Democraticness", "State A Capabilities", "State B Capabilities",
                               "Relative Capabilities", "Alliance Score"),
          add.lines=list(c("Observations", obs1, obs2, obs3, obs4),
                         c("McFadden Pseudo-R$^2$", pseudoR2.1, pseudoR2.2, pseudoR2.3, pseudoR2.4),
                         c("Cluster Robust SEs", "$\\checkmark$", "", "$\\checkmark$", ""),
                         c("Country Fixed Effects", "" , "$\\checkmark$", "" , "$\\checkmark$")),
          dep.var.labels = c("Failure to achieve goals"),
          title = "Adversary Initiated Crises",
          notes.append=FALSE, no.space=T, single.row = T,
          notes = c("Robust standard errors are clustered by crisis.$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}"))

#########################################################################
##Multinomial Logit Model################################################
#########################################################################

rm(list = ls())
load("tseries.RData")

d$crisis.mn <- "nocrisis"
d$crisis.mn[d$crisis_initA==1] <- "crisis_success"
d$crisis.mn[d$crisis_init_failA==1] <- "crisis_failure"
d$crisis.mn <- as.factor(d$crisis.mn)
d$crisis.mn <- relevel(d$crisis.mn, ref = "crisis_success")
d$crisis.mn <- relevel(d$crisis.mn, ref = "nocrisis")

var <- c("crisis.mn", "statea", "type", "age", "gender", "milservice", "combat", "hos_experience", "gwf_military", "gwf_personal", "regtransA", "polity2A", "polity2B", "LcincA", "LcincB", "Lcinc_rel", "s_wt_atop", "time", "time2", "time3")
d2 <- subset(d, select=var)
d2$choiceid <- 1:nrow(d2)
d2a <- dfidx(d2, shape = "wide", sep = "_",
             choice = "crisis.mn", idx = list(c("choiceid", "statea")),
             idnames = c(NA, "alt"))

m_mlogit <- mlogit(crisis.mn ~ 0 | type + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + Lcinc_rel + s_wt_atop + time + time2 + time3, d2a, reflevel = "nocrisis")
stargazer(m_mlogit, digits=2, intercept.bottom=T, no.space=T, single.row=T, star.cutoffs=c(0.1,0.05,0.01))

#########################################################################
##Crisis Initiation #####################################################
#########################################################################

R2logit <- function(y,model){
  R2<- 1-(model$deviance/model$null.deviance)
  return(R2)
}

crse2 <- function(x, df){
  d.cse <- na.omit(df[, c("statea", all.vars(formula(x)))])
  fc <- d.cse$statea
  m <- length(unique(fc))
  k <- length(coef(x))
  u <- estfun(x)
  u.clust <- matrix(NA, nrow=m, ncol=k) 
  for(j in 1:k){
    u.clust[,j] <- tapply(u[,j], fc, sum)
  }
  bread <-vcov(x)
  c.vcov <- bread %*% ((m / (m-1)) * t(u.clust) %*% (u.clust)) %*% bread
  result <<- coeftest(x, c.vcov)
  return(result)
}

#Model 1
m1 <- glm(crisis_initA ~ type + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + Lcinc_rel + s_wt_atop + time + time2 + time3, 
          family=binomial(link = "logit"), data = d)
pseudoR2.1 <- round(R2logit(crisis_initA, m1),3)
obs1 <- nobs(m1)
m1 <- crse2(m1, d)

#Model 4
m2 <- glm(crisis_initA ~ type + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + Lcinc_rel + s_wt_atop + as.factor(statea) + time + time2 + time3, 
          family=binomial(link = "logit"), data = d)
pseudoR2.2 <- round(R2logit(crisis_initA, m2),3)
obs2 <- nobs(m2)
m2 <- crse2(m2, d)

stargazer(m1, m2, type = "latex", omit=c("statea"),
          omit.stat=c("f", "ser", "ll", "aic"),
          label=c("table: initiation_model"),
          covariate.labels = c("Siloed", "Dictatorial", "Fragmented", "Leader Age", "Male", "Military Experience", "Combat Experience", "Time in Office", "Military Regime",
                               "Personalist Regime", "Regime Transition", "State A Democraticness", "State B Democraticness", "State A Capabilities", "State B Capabilities",
                               "Relative Capabilities", "Alliance Score", "Time Since Last Crisis", "Time Since Last Crisis$^2$", "Time Since Last Crisis$^3$"),
          add.lines=list(c("Observations", obs1, obs2),
                         c("McFadden Pseudo-R$^2$", pseudoR2.1, pseudoR2.2),
                         c("Cluster Robust SEs", "$\\checkmark$", ""),
                         c("Country Fixed Effects", "" , "$\\checkmark$")),
          dep.var.labels = c("crisis initiation"),
          title = "Crisis Initiation",
          notes.append=FALSE, no.space=T, single.row = T,
          notes = c(""))

#########################################################################
##Original ICB Values####################################################
#########################################################################

rm(list = ls())
load("analysis_original.RData")

R2logit <- function(y,model){
  R2<- 1-(model$deviance/model$null.deviance)
  return(R2)
}

crse <- function(x, df){
  d.cse <- na.omit(df[, c("crisno", all.vars(formula(x)))])
  fc <- d.cse$crisno
  m <- length(unique(fc))
  k <- length(coef(x))
  u <- estfun(x)
  u.clust <- matrix(NA, nrow=m, ncol=k) 
  for(j in 1:k){
    u.clust[,j] <- tapply(u[,j], fc, sum)
  }
  bread <-vcov(x)
  c.vcov <- bread %*% ((m / (m-1)) * t(u.clust) %*% (u.clust)) %*% bread
  result <<- coeftest(x, c.vcov)
  return(result)
}

##Model 1: 
m1.1 <- glm(crisis_init_failA ~ type, 
            family=binomial(link = "logit"), data = d.original[d.original$crisis_initA==1,])
pseudoR1.1 <- round(R2logit(crisis_init_failA, m1.1),3)
obs1.1 <- nobs(m1.1)
m1.1 <- crse(m1.1, d.original[d.original$crisis_initA==1,])

##Model 2: 
m1.2 <- glm(crisis_init_failA ~ type + as.factor(statea), 
            family=binomial(link = "logit"), data = d.original[d.original$crisis_initA==1,])
pseudoR1.2 <- round(R2logit(crisis_init_failA, m1.2),3)
obs1.2 <- nobs(m1.2)
m1.2 <- crse(m1.2, d.original[d.original$crisis_initA==1,])

##Model 3:
m1.3 <- glm(crisis_init_failA ~ type + age +gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop, 
            family=binomial(link = "logit"), data = d.original[d.original$crisis_initA==1,])
pseudoR1.3 <- round(R2logit(crisis_init_failA, m1.3),3)
obs1.3 <- nobs(m1.3)
m1.3 <- crse(m1.3, d.original[d.original$crisis_initA==1,])

##Model 4:
m1.4 <- glm(crisis_init_failA ~ type + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop + as.factor(statea), 
            family=binomial(link = "logit"), data = d.original[d.original$crisis_initA==1,])
pseudoR1.4 <- round(R2logit(crisis_init_failA, m1.4),3)
obs1.4 <- nobs(m1.4)
m1.4 <- crse(m1.4, d.original[d.original$crisis_initA==1,])

stargazer(m1.1, m1.2, m1.3, m1.4, type = "latex", omit=c("statea"),
          omit.stat=c("f", "ser", "ll", "aic"),
          label=c("table: main_model"),
          covariate.labels = c("Siloed", "Dictatorial", "Fragmented", "Leader Age", "Male", "Military Experience", "Combat Experience", "Time in Office", "Military Regime",
                               "Personalist Regime", "Regime Transition", "State A Democraticness", "State B Democraticness", "State A Capabilities", "State B Capabilities",
                               "Relative Capabilities", "Alliance Score"),
          add.lines=list(c("Observations", obs1.1, obs1.2, obs1.3, obs1.4),
                         c("McFadden Pseudo-R$^2$", pseudoR1.1, pseudoR1.2, pseudoR1.3, pseudoR1.4),
                         c("Cluster Robust SEs", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                         c("Country Fixed Effects", "" , "$\\checkmark$", "" , "$\\checkmark$")),
          dep.var.labels = c("Failure to Achieve Goals"),
          title = "National Security Institutions and International Crisis Performance",
          notes.append=FALSE, no.space=T, single.row = T,
          notes = c("Robust standard errors are clustered by crisis.$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}"))
